We need doxygen comments for all descriptions of non-override methods and for all member variables. Three slash comments for the variables, long form comments for methods, including @param and @return tags.

It's going to be important to implement this for anything off-axis.

By this I mean that if the bar is oriented in a direction off the Cartesian axes, this will be nonzero. You can ignore this for the first cut, but this will be very important to get right.

Should be set = 0?

I guess that is what this accomplishes. It just looks strange. I'd replace it with an =

Lines 156 and 271 don't seem consistent.

What would be a better way to get the element length/volume?

These are all computed assuming that there is a penalty stiffness in the axial and perpendicular directions, which isn't the case anymore with the bond slip model. You would need to apply the penalty stiffness in the perpendicular directions, but subtract that off in the axial direction and replace it with the tangent from the bond slip model. We do that kind of thing in contact in the MechanicalContactConstraint class. Same idea here. I imagine this is the main source of the convergence trouble.

I was just thinking -- there's a newly introduced Constraint in tensor_mechanics that Swetha wrote that does something very similar to this Constraint: NodalFrictionalConstraint. It represents a frictional slider that uses a penalty constraint until it slips. I just looked at the code for the Jacobian in that class. It computes the nonlinear Jacobian, but it doesn't have any of the code I was expecting to see to deal with the directionality of the Jacobian. That's because the NodalFrictionalConstraint class is only applied to the dof in the axial direction of the beam, and it assumes that the beam is aligned in one of the Cartesian directions. Another separate Constraint (NodalStickConstraint) deals with tying it in the other two directions. Your constraint is dealing with all directions (which is probably preferable).

The bondslip contribution towards the residual is constant. So I only considered the penalty contributions towards jacobian calculation. I believe I am using the normal direction to assign the contribution. If the truss member is inclined, only then the penalty would apply to both x- and y- directions.